// CG CO2 dimer model will use NVT-MC to parameterize LJ against bulk radial distribution
//This version uses NPT ensemble sampling
//revised: 6.23.2011
//S.Stauffer (writer) 

#include <iostream> 
#include <fstream> 
#include <string> 

using namespace std; 

#include "mt.c" 
#include "CO2_constants.h"
#include "CO2_Initialize.h"	// INITIALIZATION CODE 
#include "FileIO.cc"
#include "Configuration.cc" 
#include "CO2_NVT_kernel.cu" 	// GPU KERNEL CODE 
#include "CO2_NPT.h"		// NPT 
//#include "CO2_Virial.h"		// VIRIAL 

void monteCarlo(int nSteps) {

	int ref = 1;
	int  moveCountNVT=0, acceptCountNVT=0;

	prec x0w;  //e_kin, pot2,acceptRatiol,x0wl,fmw;
	
	srand(time(NULL));

	nvtLog.open("NVT-CO2.dat"); 
	nvtLog << "step \t" << "energy (kCal) \t" << "move size (A) \t" << "acceptance ratio" << endl; 

	nptLog.open("NPT-CO2.dat");
	nptLog << "step \t" << "density \t" << "average density \t" << "energy (kCal) \t" << "acceptance ratio \t" << "box length \t" << "enthalpy change" << endl;

	for(int stepCount=0;stepCount<nSteps;stepCount++) {

		//Random number generation is selecting all particles equally
		moveSelect = ((int)(genrand_real()*nMolecules ));  

		if (moveSelect < 0) moveSelect = 0;

		moveCountNVT++; 
		
		ref = 1;	

		savePositions("save"); 

		prec initPotential = calculatePotential(ref, moveSelect, hAtom); 

		translatePositions(); 			

		prec propPotential = calculatePotential(ref, moveSelect, hAtom);
		prec deltaPotential = propPotential - initPotential; 
	   
		prec r = genrand_real();

		if ((exp(-1*deltaPotential/(BOLTZMANN_C*tempK)))>r) {

			acceptCountNVT++;
			potentialEnergy += deltaPotential;
		} else { 
			savePositions("restore"); 
		} 

		prec acceptRatio=(acceptCountNVT)/(double)(moveCountNVT);

		if (simulationType == "EQ" && (stepCount+1)%2000==0) {

			x0w=(acceptRatio-0.4)/0.4;
			moveSize=moveSize+x0w*moveSize;
			dRm = dRm+ x0w*dRm;

			if (stepCount == (nSteps-1)) { 
				writeConfiguration(); 
			} 
		}

        if (simulationType == "PROD" && stepCount%2500==0) {

			writeConfiguration(); 
		} 

		if (stepCount % 10000 == 0) { 

			nvtLog << stepCount << "\t\t" << potentialEnergy/ENERGY_CF << "\t\t" << moveSize << "\t\t" << acceptRatio << endl;
		}

		if (nptSelect == "ON" && stepCount%250 == 0) NPT(); 

		cout << stepCount << "\t" << potentialEnergy/ENERGY_CF << endl; 

	
	} // End of MC iterations loop
} 


int main(int argc, char *argv[]) {

    cudaSetDevice(deviceID);					// Initialize the GPU

	long s=-12356;								// Initialize the RNG
	init_genrand(s);

    initSimulation();							// Initialize simulation parameters
												// Perform initialization BEFORE allocation in order to properly set nMolecules for dynamic allocation
    allocateMemory();							// Allocate memory on host and device
	
	readConfiguration();						// Read configuration file
	
	normalizePositions();						// Place everything into the box

	potentialEnergy = calculatePotential(0, 0, hAtom); // Calculate reference potential

	cout << "Reference Energy (J):" << potentialEnergy/ENERGY_CF <<"\n";

	monteCarlo(nSteps);

	freeMemory();								// Free memory on host and device
}

